#include "sphere.h"
#include "../core/transform.h"
#include "../core/sampling.h"
#include "../core/geometry.h"

Sphere::Sphere(const Transform *o2w, const Transform *w2o, bool ro, Reference<Material> mat, float rad) : 
	Shape(o2w, w2o, ro, mat), radius(rad) {
}

BBox Sphere::ObjectBound() const {
	return BBox(Point3(-radius, -radius, -radius), Point3(radius, radius, radius));
}

bool Sphere::IntersectP(const Ray &r) const {
	// Transform Ray to object space
	Ray ray;
	(*WorldToObject)(r, &ray);
	// Compute quadratic sphere coefficients
	float A = ray.d.LengthSqured();
	float B = 2 * Dot(ray.d, (Vector3)ray.o);
	float C = Vector3(ray.o).LengthSqured() - radius * radius;
	// Solve quadratic equation
	float t1, t0;
	if(!Quadratic(A, B, C, &t0, &t1))
		return false;

	if(t0 > ray.maxt || t1 < ray.mint)
		return false;
	if(t0 < ray.mint && t1 > ray.maxt)
		return false;
	
	return true;
}

bool Sphere::Intersect(const Ray &r, Intersection *isect) const {
	// Transform Ray to object space
	Ray ray;
	(*WorldToObject)(r, &ray);
	// Compute quadratic sphere coefficients
	float A = ray.d.LengthSqured();
	float B = 2 * Dot(ray.d, (Vector3)ray.o);
	float C = Vector3(ray.o).LengthSqured() - radius * radius;
	// Solve quadratic equation
	float t1, t0;
	if(!Quadratic(A, B, C, &t0, &t1))
		return false;
	if(t0 > ray.maxt || t1 < ray.mint)
		return false;
	float thit = t0;
	if(t0 < ray.mint) {
		thit = t1;
		if(t1 > ray.maxt)
			return false;
	}
	// Update r.maxt
	r.maxt = thit;
	isect->t = thit;
	isect->shape = this;
	return true;
}

void Sphere::FillIntersection(const Ray &ray, Intersection *isect) const {
	isect->p = ray(isect->t);
	Vector3 local = (*WorldToObject)((Vector3)isect->p - Vector3(0.0, 0.0, 0.0));
	float theta = acosf(Clamp(local.z/radius, -1.f, 1.f));
	float phi = atan2f(local.y, local.x);

	if (phi < 0)
		phi += 2*M_PI;

	isect->uv.x = phi * INV_TWOPI;
	isect->uv.y = theta * INV_PI;

	Vector3 dpdu, dpdv;
	dpdu = Vector3(-local.y, local.x, 0) * (2*M_PI);
	isect->dpdu = (*ObjectToWorld)(dpdu);
	isect->geometryFrame.n = Normalize((*ObjectToWorld)((Normal)local));
	float zrad = std::sqrt(local.x*local.x + local.y*local.y);
	float cosPhi, sinPhi;
	if (zrad > 0) {
		float invZRad = 1.0f / zrad;
		cosPhi = local.x * invZRad,
		sinPhi = local.y * invZRad;
		dpdv = Vector3(local.z * cosPhi, local.z * sinPhi, -sinf(theta)*radius) * M_PI;
		isect->dpdv = (*ObjectToWorld)(dpdv);
		isect->geometryFrame.s = Normalize(isect->dpdu);
		isect->geometryFrame.t = Normalize(isect->dpdv);
	} else {
		// avoid a singularity
		cosPhi = 0, sinPhi = 1;
		dpdv = Vector3(local.z * cosPhi, local.z * sinPhi, -sin(theta)*radius) * M_PI;
		isect->dpdv = (*ObjectToWorld)(dpdv);
		CoordinateSystem(isect->geometryFrame.n, &isect->geometryFrame.s, &isect->geometryFrame.t);
	}

	// Compute dndu dndv

	Vector3 d2Pduu = - 4 * M_PI * M_PI * Vector3(local.x, local.y, 0.);
	Vector3 d2pduv = 2 * M_PI * M_PI * local.z * Vector3(-sinPhi, cosPhi, 0.);
	Vector3 d2pdvv = - M_PI * M_PI * Vector3(local.x, local.y, local.z);

	float E = Dot(dpdu, dpdu);
	float F = Dot(dpdu, dpdv);
	float G = Dot(dpdv, dpdv);
	Vector3 N = Cross(dpdu, dpdv);
	float e = Dot(N, d2Pduu);
	float f = Dot(N, d2pduv);
	float g = Dot(N, d2pdvv);

	float invEGF2 = 1.f / (E*G - F*F);
	Normal dndu = (Normal)((f*F - e*G) * invEGF2 * dpdu +
					(e*F - f*E) * invEGF2 * dpdv);
	Normal dndv = (Normal)((g*F - f*G) * invEGF2 * dpdu +
					(f*F - g*E) * invEGF2 * dpdv);

	isect->dndu = (*ObjectToWorld)(dndu);
	isect->dndv = (*ObjectToWorld)(dndv);

	if(ReverseOrientation)
		isect->geometryFrame.n *= -1.f;

	isect->shadingFrame = isect->geometryFrame;

	isect->shape = this;
	isect->hasUVPartials = false;

	isect->objectToWorld = this->ObjectToWorld;
	isect->worldToObject = this->WorldToObject;

	isect->isLight = (this->arealight != NULL);
	if(isect->isLight) {
		//isect->le = this->arealight->GetLe();
		Warning("isect->le not initialized yet\n");
	}
}
float Sphere::Area() const {
	return 4.f * M_PI * radius * radius;
}
Point3 Sphere::Sample(const Point2 &sample, Normal *Ns) const {
	Point3 p = Point3(0,0,0) + radius * UniformSamplingSphere(sample);
	*Ns = Normalize((*ObjectToWorld)(Normal(p.x, p.y, p.z)));
	if(ReverseOrientation) 
		*Ns *= -1.f;
	return (*ObjectToWorld)(p);
}

/*
* Improved sampling strategy given in
* "Monte Carlo techniques for direct lighting calculations" by
* Shirley, P. and Wang, C. and Zimmerman, K. (TOG 1996)
*/
Point3 Sphere::Sample(const Point3 &p, const Point2 &sample, Normal *Ns) const {
	Point3 Pcenter = (*ObjectToWorld)(Point3(0.,0.,0.));
	if(DistanceSqured(p,Pcenter) - radius * radius < EPSILON(4)) {
		// P is inside the sphere
		// Uniform sampling
		return Sample(sample, Ns);
	}
	Vector3 wc = Normalize(Pcenter - p);
	Vector3 wcX, wcY;
	CoordinateSystem(wc, &wcX, &wcY);
	float cosThetaMax = sqrtf(max(0.f, 1.f - radius*radius/DistanceSqured(Pcenter, p)));
	Intersection isect;
	Point3 ps;
	Ray r(p, UniformSamplingCone(sample, cosThetaMax, wc, wcX, wcY));
	if(!Intersect(r, &isect))
		ps = Pcenter - radius * wc;
	else
		ps = r(isect.t);
	*Ns = (Normal)Normalize(ps - Pcenter);
	if(ReverseOrientation)
		*Ns *= -1.f;
	return ps;
}

float Sphere::Pdf(const Point3 &p, const Vector3 &wi) const {
	Point3 Pcenter = (*ObjectToWorld)(Point3(0.,0.,0.));
	if(DistanceSqured(p,Pcenter) - radius * radius < EPSILON(4)) {
		// P is inside the sphere
		// Uniform sampling
		return Shape::Pdf(p, wi);
	}
	float cosThetaMax = sqrtf(max(0.f, 1.f - radius*radius/DistanceSqured(Pcenter, p)));
	return UniformSamplingConePdf(cosThetaMax);
}
